Energy landscape analysis and time-series clustering analysis of patient state multistability related to rheumatoid arthritis drug treatment: The KURAMA cohort study

Rheumatoid arthritis causes joint inflammation due to immune abnormalities, resulting in joint pain and swelling. In recent years, there have been considerable advancements in the treatment of this disease. However, only approximately 60% of patients achieve remission. Patients with multifactorial diseases shift between states from day to day. Patients may remain in a good or poor state with few or no transitions, or they may switch between states frequently. The visualization of time-dependent state transitions, based on the evaluation axis of stable/unstable states, may provide useful information for achieving rheumatoid arthritis treatment goals. Energy landscape analysis can be used to quantitatively determine the stability/instability of each state in terms of energy. Time-series clustering is another method used to classify transitions into different groups to identify potential patterns within a time-series dataset. The objective of this study was to utilize energy landscape analysis and time-series clustering to evaluate multidimensional time-series data in terms of multistability. We profiled each patient’s state transitions during treatment using energy landscape analysis and time-series clustering. Energy landscape analysis divided state transitions into two patterns: “good stability leading to remission” and “poor stability leading to treatment dead-end.” The number of patients whose disease status improved increased markedly until approximately 6 months after treatment initiation and then plateaued after 1 year. Time-series clustering grouped patients into three clusters: “toward good stability,” “toward poor stability,” and “unstable.” Patients in the “unstable” cluster are considered to have clinical courses that are difficult to predict; therefore, these patients should be treated with more care. Early disease detection and treatment initiation are important. The evaluation of state multistability enables us to understand a patient’s current state in the context of overall state transitions related to rheumatoid arthritis drug treatment and to predict future state transitions.


Introduction
Rheumatoid arthritis (RA) causes joint inflammation due to immune abnormalities, resulting in joint pain and swelling.In recent years, there have been considerable advancements in the treatment of RA, partly due to the development of drugs such as methotrexate (MTX), biologic disease-modifying anti-rheumatic drugs (bDMARDs), and targeted synthetic DMARDs (tsDMARDs) such as Janus kinase (JAK) inhibitors; furthermore, the "treat-to-target (T2T) algorithm", in which treatment is periodically adjusted to a target disease index, has led to improvements in RA treatment [1][2][3][4][5].However, even with these approaches, only approximately 60% of patients achieve remission.Therefore, 10-20% of RA patients who are treatment-refractory have been identified as having difficult-to-treat (D2T) RA.In 2020, the European League Against Rheumatism (EULAR) published the D2T RA EULAR Definition.The development of appropriate treatments for refractory patients is urgently needed [6][7][8].
Due to advances in information technology, a variety of digitized data from daily practice, such as electronic medical records and various laboratory test values, can be collected.Therefore, expectations regarding "real-world evidence" are increasing [9][10][11].However, the aim of many clinical studies is to consider the whole treatment as a "single intervention" and to determine the effectiveness and safety of that intervention while eliminating bias as much as possible.Multidimensional time-series data over the entire course of treatment are rarely analyzed in a temporal manner.In contrast, daily practice requires that the patient's condition be observed over time and that the treatment method be selected in a timely manner to optimize efficacy and survival time.In other words, in daily practice, the whole treatment is not considered a "single intervention."Instead, the optimal treatment is selected based on the constantly changing state of the patient.To obtain high-quality real-world evidence, it is necessary to develop a method for analyzing multidimensional time-series data in a time-dependent manner over the course of treatment and providing information that is more consistent with decision-making in daily practice [12,13].
In multidimensional time-series data analysis in medicine, there are examples of applying dynamic treatment regimens (DTRs) and deep learning.However, establishing DTRs requires data on all treatments beginning at the initial visit in addition to an intermediate covariate history.The computational requirements for model building are high.Therefore, DTRs are often used in clinical trials or observational studies where multiple treatment regimens that are predicted in advance are compared in a time-dependent manner [14][15][16][17].Furthermore, the logic of deep learning is not very transparent.In many cases, the rationale is unclear, thus making it difficult to implement deep learning in daily practice [18][19][20][21][22][23][24].
The Kyoto University Rheumatoid Arthritis Management Alliance (KURAMA) cohort database of the Rheumatology Center of Kyoto University Hospital includes data on initial consultation, follow-up, blood tests, etc., from all patients with RA visiting the center.The database has both clinical and research applications.It includes information on more than 3,000 patient registrations with more than 40,000 disease activity data points and has already been used extensively for the analysis of drug administration and disease activity [25][26][27][28].The KURAMA cohort database includes high-quality, multidimensional time-series data and is considered to be suitable for time-series data analysis.
Energy landscape analysis is used to estimate a model distribution of maximum likelihood from data, making it possible to visualize the ease with which transitions occur between highand low-energy states and thus enabling intuitive interpretation [29,30].Patients with chronic or multifactorial diseases shift between states from day to day.These patients may remain in a good or poor state with little change, or their condition may vary frequently.Energy landscape analysis can quantitatively evaluate the stability/instability of each state as energy, which is not possible with the conventional method of simply clustering state variables.By assessing stable/ unstable states in addition to conventional good/poor states, energy landscape analysis makes it possible to determine when interventions are effective and which states are difficult to treat.Energy landscape analysis is often used in areas such as protein folding and stability analysis [31][32][33] and is considered useful for visualizing the state transitions of patients in real-world practice.
Time-series clustering is another method used to classify transitions into different groups to identify potential patterns within a time-series dataset [34][35][36][37][38][39].Dynamic time warping (DTW) is a type of time-series clustering that measures the distance and similarity between time-series data by finding the shortest path, which can be identified by summing the distance (i.e., the absolute value of the error) between the points of two time series.DTW can determine the similarity of time series even if the length and period of the time series are different [35,40].
We hypothesized that using energy landscape analysis and time-series clustering with DTW to evaluate the multidimensional time-series data in the KURAMA cohort would enable the visualization of transitions between time-dependent states in patients with RA during drug treatment, thus providing useful information for achieving the goals of RA treatment.
The purpose of this study was to utilize energy landscape analysis and time-series clustering with DTW to evaluate multidimensional time-series data in the KURAMA cohort in terms of multistability, thereby facilitating the achievement of RA treatment goals.

Study design and participants
This single-center, retrospective, observational study utilized the KURAMA cohort database (S1 Checklist).The study protocol adhered to the principles of the Declaration of Helsinki and was approved by the Kyoto University Graduate School of Medicine Medical Ethics Committee through a central collective review (R2820), and written informed consent for study participation was obtained from all patients.The study received approval on March 30, 2021, granting access to the data from that point onward.
All patients who presented to the Rheumatology Center of Kyoto University Hospital and who met the 1987 or 2010 RA classification criteria [41] were enrolled in the KURAMA cohort study, and clinical and functional data were recorded at baseline and at each visit during the study.The inclusion criteria were as follows: patients with RA enrolled between January 1, 2011, and December 31, 2018; no previous medication history at the first visit; and onset of clinical remission within 3 years or follow-up for up to 3 years.

Variables
We defined a model of patient state transitions in RA drug treatment based on T2T.The patient's state shifts from time to time.Fig 1 shows the ease and direction of transition of the patient's state based on high and low energy in RA practice.
When energy is low, the patient is stable in a good or poor state.On the other hand, when the energy is high, the patient is unstable and transitions easily to another state.The highenergy state is considered to indicate the period of effective treatment.Individual states were evaluated using a comprehensive disease activity assessment."Good stability" was defined as meeting functional remission criteria based on the Health Assessment Questionnaire (HAQ), and "poor stability" was defined as falling below the energy threshold without functional remission (i.e., treatment dead-end) [42].The term "treatment dead-end" indicates that the patient will be not in remission, regardless of all potential future treatment sequences.In this study, the treatment period was 3 years, and the state of each patient was evaluated up to 12 times.In RA practice, the physician collects information about the patient's state through visual examination, palpation, blood tests, and imaging tests.Variables related to disease activity, bone destruction, and immunological response are used in comprehensive disease activity assessments in daily practice and were used to characterize the patient's state in this study, including: • Rheumatoid factor (RF)  [43] • Swollen joint count-28 joints (SJC28) • Tender joint count-28 joints (TJC28) The time-series data that served as input values for the energy landscape analysis were binarized as high activity (i.e., nonremission) = 1 and low activity (i.e., remission) = -1.RF and ESR1h were binarized based on blood test reference values.PtVAS, DrVAS, SJC28, and TJC28 were binarized based on Boolean remission criteria [44], and STAGE was binarized based on Steinbrocker's staging classification [43].The remission criteria were as follows: • RF = 15 #Unit: IU/mL • ESR1h (male) = 10 #units: mm • ESR1h (female) = 20 #units: mm The higher the energy is, the more unstable the state and the easier it is to obtain a therapeutic effect.At initial consultation, many patients present in a poor state and with high disease activity (i.e., high energy).With treatment, they may transition to a good state or stabilize in a poor state.The physician observes the patient's state and executes a treatment strategy to stabilize the patient in a good state or to prevent stabilization in a poor state below a certain energy threshold.The objective of this study was to visualize the state transitions of these patients as a population.Since there are two patterns per factor (i.e., +1/-1), analyzing seven factors generates 128 (i.e., 2^7) activity patterns.The physician selected the optimal treatment according to this information (i.e., state number) (Fig 2).

Analytical methods
Overview of the analytical methods.Energy landscape analysis was used to assess and visualize state multistability as "energy."Dynamic systems were formulated with a Boltzmann machine that used the Boltzmann distribution with energy to define the probability distribution.The Boltzmann machine reduced multidimensional data to a measurement considered to represent "energy."An index of multistability was assigned to each state to identify whether the disease condition was improving or worsening.Furthermore, multidimensional timeseries clustering of patients enabled us to classify and visualize the time-series transitions of unstable and stable states.
Energy landscape analysis and the Boltzmann machine.For each patient (s), there are seven different factor values (i) observed at each of the 12 time points (t) every 3 months over 3 years.Let O ¼ fTJC28; SJC28; STAGE; DrVAS; PtVAS; ESR1h; RFg and T = {1,2,. ..,11,12}.Each factor x i (t) is binarized (-1 or +1) according to the clinical criteria for i2O and t2T.We apply a Boltzmann machine [45] that extends the deterministic dynamics of the Hopfield model, a well-known model of associative memory, to stochastic dynamics.
The Boltzmann machine is defined as a multidimensional Boltzmann distribution.The definition of the distribution includes energy.A stochastic model on an undirected graph G (O, E) is defined, where E is the set of (i,j) links between nodes i,j2O.
where x i = −1 or +1, vector θ = {θ i }, and matrix W = {w ij } for i,j2O.G (O, E) are assumed in the link structure of the Boltzmann machine, such that w ij = w ji and w ii = 0.The third equation is called the energy function.The second equation is the normalizing constant of the first equation.According to the definition of the first equation, this Boltzmann machine is a mathematical model in which the energy function is smaller than the monotonically increasing nature of the exponential function, and the activity pattern appears with higher probability.The parameters θ and W are trained to match the probability of occurrence of the actual observed data.For each person s2S, the observation vector obtained at each time t2T is assumed to be generated from an independent homoscedastic distribution.For the obtained dataset D ¼ fx s i ðtÞji 2 O; t 2 T; s 2 Sg, the maximum likelihood method is applied to estimate the parameters s θ and W that satisfy the following equation: The parameter that maximizes this log-likelihood function, log L D (θ,W), is obtained using the gradient ascent method as follows: where E old is the expected value from the Boltzmann distribution using the previous parameter in the dataset, X(i) is the i-th element of X,i,j2O, |S| is the number of elements in S, and ε = 0.2 is the learning rate.Up to 5,000,000 iterations were performed.Disconnectivity graph.First, a minimal activity pattern X = {x i |i2O} is defined such that for any activity pattern Y with one different node in X, F(Y)≧F(X) holds for the given parameters s θ and W. Next, a path a with X and Y is defined.Let A(X,Y,a) = {Z|Z transform X one node at a time until it can be transformed into Y}.Note that this set exists for the number of paths a. Finally, the energy of the highest hill to be surmounted by path a is max AðX;Y;aÞ FðZÞ.
Time-series clustering.We chose the variables y id (t) based on the activity pattern and chose the energy and HAQ for patient id and time t.Clustering was performed using the kmeans method with distances based on the dynamic contraction method toward time-series vectors fy id ðtÞg t2T , where y id ðtÞ 2 R jOj � R 2 ; for id = 1,2,� � �,N, T is the time point, and N is the number of samples [48].
The number of clusters K was determined by the elbow method or silhouette analysis [49,50].The initial central value fc id ðtÞg t2T of cluster k was randomly chosen for c id ðtÞ 2 R jOj � R 2 .The following method of minimizing the within-cluster sum of squared errors of prediction (SSE) was used.Iterations by the k-means method were performed up to 50 times.The center value was updated by the mean vector in the cluster.for u<S and v<S, where |・| is the Euclidean distance.The DTW algorithm [51] for calculating the distance between two time series uses least-cost elasticity matching, which does not allow the time series to intersect.

Analysis environment and preprocessing
Oracle Autonomous Data Warehouse and Oracle Cloud Infrastructure Data Science were used as analysis environments.The Statistics and Machine Learning Toolbox and the Energy Landscape Analysis Toolbox v1.2 [29] of MATLAB R2022b were used for energy landscape analysis.
For missing value completion, STAGE, which comprises categorical data, was substituted for the before and after data.All other factors were numerical data, and linear interpolation of time-series data was used; for patients with fewer than 12 points, data were generated in the same way as for missing value imputation above.For example, for patients who entered remission or who withdrew from the study, data were generated in the same way as for the final point until the patients reached 3 years of remission.In the case of missing time-series points, categorical STAGE data were assigned before and after the values, and the numerical data were linearly interpolated.

Participants and descriptive data
A flowchart of the subject selection, the 3-year trends for the seven factor values, and the distribution of the drug therapies administered are shown in Fig 3.
From the 3,439 individuals enrolled in the KURAMA cohort from 2011 to 2018, we excluded 2,401 individuals who were duplicates or who lacked laboratory values or medication information within 3 years of their first visit.In addition, we excluded 441 patients with a history of past medication use at the time of their first visit.Thus, 597 patients were ultimately included in the analysis.
The baseline characteristics of the 597 subjects are shown in Table 1.

Results of the energy landscape analysis
The parameters of the Boltzmann machine were as follows: θ = -0.4259-0.2, and S1 Table ) revealed that patient state could be divided into two fixed disease state patterns: "good stability leading to remission" (blue) and "poor stability (i.e., treatment dead-end) (red).The higher the energy value is, the more unstable the state is and the easier it is to transition to another state.On the other hand, the lower the energy value is, the more fixed the state is and the harder it is to transition to another state.In other words, higher energy values indicate that the patient is more responsive to drug treatment or that the rate of disease deterioration exceeds the effect of drug treatment.On the other hand, lower energy levels mean that the patient is more likely to remain in the current state and that the disease is less likely to improve or worsen.
Fig 4B shows whether each of the seven factors met the remission and nonremission criteria according to the state numbers that take the minimal energy of the patterns of good stability and poor stability.The factors shown in black met the remission criteria, while those in white did not.In the "good stability" pattern, all factors except RF and PtVAS met the remission criteria.On the other hand, in the "poor stability" pattern, none of the factors reached the remission criteria.
The patient's state could switch between the two patterns as a result of drug treatment effects or other factors.The quadrant in which the energy is lower than the threshold and the activity pattern is "good stability" is G-L; the quadrant in which the energy is higher than the threshold and the activity pattern is "good stability" is G-H; the quadrant in which the energy is higher than the threshold and the activity pattern is "poor stability" is P-H; and the quadrant in which the energy is lower than the threshold and the activity pattern is "poor stability" is P-L.
Fig 5B shows the number of patients in each quadrant, with a marked increase in G-L and a decrease in P-L until approximately two points (i.e., 6 months), followed by a gradual increase or decrease.P-L decreased until approximately two points (i.e., 6 months), with no significant change thereafter.G-H decreased significantly throughout the entire period.G-H did not increase or decrease significantly across the entire treatment period.
Fig 5C shows the HAQ scores in each quadrant: G-L patients met the remission criteria for almost the entire period; the number of P-L patients who met the remission criteria decreased over time; the number of G-H and P-H patients who met the remission criteria changed     minimally over time; and the number of P-H patients who met the criteria did not change over time.
In terms of overall treatment, the number of patients who responded to drug therapy and who demonstrated an improved disease status increased markedly until approximately 6 months after the start of treatment, after which the rate of increase gradually slowed.The rate of change plateaued after 1 year.

Results of time-series clustering
Individual patient state transitions are highly variable, making it difficult to identify patterns across the population (S1 Appendix).Time-series clustering was performed based on activity  patterns, energy, and HAQ scores, and patient profiling was performed for each cluster.The number of clusters was set to "3" based on comprehensive analysis using the elbow method and silhouette technique [48,49] (Fig 6 and Table 3).Cluster 0 (toward good stability) (blue) was defined as patients with low energy, 80% of whom eventually assumed a good stability pattern and became stable.Cluster 1 (toward poor stability) (red) was defined as patients with low energy, almost all of whom were characterized by the "poor stability" pattern and had high HAQ scores.Cluster 2 (unstable) (green) consisted of patients with relatively high energy, almost all of whom demonstrated a "poor stability" pattern in the early stage of treatment.For the patients in Cluster 2 (unstable)", the treatment response was not yet determined, and the majority of them met the remission criterion based on the HAQ score, although a high proportion of these patients exhibited a "poor stability" pattern in the early stage of treatment.Although the age at disease onset did not differ significantly between clusters, the age at initial diagnosis, duration of disease, and laboratory tests at initial diagnosis did differ, with the "poor stability" cluster (Cluster 1) showing a greater tendency for age and all laboratory values at initial examination than the other clusters.
The state transitions of individuals from each cluster are shown in the S2 Appendix.

Key results
The whole treatment course of patients in the KURAMA cohort was analyzed in a time-dependent manner.Energy landscape analysis revealed that state transitions were divided into two patterns: "good stability leading to remission" and "poor stability below the energy threshold without functional remission (i.e., treatment dead-end)."The energy threshold cutoff for switching between patterns was approximately -1.48.The states were aggregated into four quadrants based on this energy threshold, and the number of patients in each quadrant was calculated.
The number of patients who demonstrated status improvement increased markedly until approximately 6 months after the start of treatment, after which the rate of increase gradually slowed.The rate of change plateaued after 1 year.Patient profiling by time-series clustering showed that patients were classified into three clusters: "toward good stability," "toward poor stability," and "unstable."Patients in the "unstable" cluster are considered to have clinical courses that are difficult to predict in RA practice; therefore, these patients should be treated with more care.Age at initial diagnosis, duration of disease, and laboratory values at initial diagnosis tended to differ by cluster, with the "toward poor stability" cluster exhibiting a higher age and higher laboratory values at initial examination than the other clusters.Compared to the "toward good stability" cluster, the "unstable" cluster had a higher age at initial examination and a longer duration of illness, but the "unstable" cluster tended to have better laboratory values at initial examination.The difference between these two clusters was marginal; however, there was a possibility that the initiation of effective treatment was delayed in the "unstable" cluster due to the relatively favorable initial conditions at the time of the first visit.Further research is warranted; however, the results of this study indicate the importance of early disease detection and treatment initiation, especially in "unstable" cases.
Energy landscape analysis was performed to visualize the state transitions of the patient population as a whole.However, this approach alone did not allow for the visualization of individual patient state transitions.We utilized time-series clustering to clarify the specific state transitions occurring within the given population.These methodologies complemented each other, revealing insights into both the collective state transitions of the population and the characteristic attributes of patients within that population.Traditionally, treatment decisions have been made based on past treatment histories and the latest evaluation of disease activity.Energy landscape analysis and time-series clustering analysis could yield insights regarding how a patient's state might transition in the future and the optimal timing for effective treatment.

Limitations
This was a single-center observational study at Kyoto University Hospital that used registry data from 2011 to 2018.This study also did not evaluate the effectiveness of state transition visualization in daily RA drug therapy.

Interpretation
Currently, RA treatment practices are based on practice guidelines, and treatment approaches are determined based on evaluations of the efficacies of previously administered agents.Specific treatment strategies for individual patients are determined based on data obtained from comprehensive disease activity assessments and on the treatment history from the initial diagnosis to the present, and future state transitions rely on speculation based on physicians' clinical experience.In our study, patient response to drug therapy and improved disease status significantly increased for approximately 6 months after the start of treatment and then plateaued after 1 year.Furthermore, our findings demonstrated that the patient population could be divided into three clusters: "toward good stability," "toward poor stability," and "unstable."We quantitatively clarified the timing and duration of significant treatment effects through time-dependent analysis of the state of patients with RA, which is ever-changing due to drug therapy.In addition to the conventional assessment of overall disease activity, the evaluation of energy (i.e., state multistability) enables us to understand the relevance of the current state in the overall state transition related to RA drug treatment and to predict future state transitions.In the realm of artificial intelligence (AI) applications in healthcare, there is a fascinating parallel with the domain of chess, shogi, and go AI, where professionals occasionally find themselves defeated by seemingly unorthodox moves beyond human comprehension [52].Similarly, in the context of real-world practice visualization derived from our methodologies, the recommended treatment may significantly deviate from established clinical guidelines.Despite low disease activity, AI could suggest early and intensive administration of biologic agents.In such cases, it becomes imperative to assume that the proposed therapeutic approach of AI holds merit and to therefore perform rigorous validation through randomized controlled trials (RCTs).Our research suggests the emergence of an innovative approach to the future landscape of medical research and practice, tentatively termed "AI-based RCTs," to systematically investigate and corroborate the efficacy of AI-suggested clinical interventions.In the context of our study, for patients who are likely to experience treatment dead-end, it may be necessary to select a more effective treatment, such as a biologic, at an earlier stage.For patients whose condition improves, it may be possible to offer drug discontinuation or dose reduction.Clinical research involving "unstable" patients may allow for more effective clinical evaluation with smaller numbers of patients.Furthermore, state transition visualization may be a useful tool for facilitating communication between physicians and patients.Ultimately, by selecting treatments according to a patient's cluster and their state transitions, it may be possible to adopt a personalized medicine approach targeting each patient's state.This study suggested the possibility of optimizing the treatment plan based on the whole treatment course.

Generalizability
This was a single-center observational study.Kyoto University Hospital is a tertiary hospital and may accommodate patients with more severe illness than general hospitals and clinics; therefore, a multicenter study is needed.In addition, since this study was based on data collected in 2018, new drugs such as JAK inhibitors were likely to be used less frequently.Therefore, data from 2019 and beyond should be examined in future studies.Analyzing more recent data will enable the visualization of state transitions of RA drug therapy that more accurately reflect the reality of daily practice.
We believe that our proposed approach is useful as a visualization method for multidimensional time-series data in medicine and can be applied to diseases other than RA.In general, it is difficult to intuitively interpret multidimensional time-series data in medicine.We have made it possible to readily interpret state transitions in the drug treatment of patients with RA by consolidating the seven factors used in comprehensive disease activity assessment in daily practice into two dimensions, namely, state number and energy, and clustering along the time course.In many diseases, even the best treatment does not lead to a cure or remission; therefore, conservative treatment is used to prevent deterioration and maintain a good state.For multifactorial diseases such as diabetes and hypertension, target treatment values such as HbA1C and blood pressure levels are set by guidelines, and treatment strategies are chosen based on concepts equivalent to T2T in RA [53,54].Various factors, such as diet, lack of exercise, and stress, affect the condition of patients with multifactorial diseases.Nonetheless, advances in digital health technology have led to increased adoption of personal health records and health-related apps; therefore, the collection of frequently measured time-series sensor data from wearable devices is becoming increasingly feasible [55][56][57].We believe that our method may provide useful information for optimizing treatment plans to achieve behavioral change and personalized medicine by reducing the dimensionality of multiple factors and visualizing state transitions.
Considering the application of this research to personalized medicine in daily practice, it is necessary to develop an application that can calculate energy based on these seven factors and display its position in the state transition pattern in real time ( Fig 7).
In addition, standardization of medical practice is necessary for high-quality multidimensional time-series data collection.Conventional medical practice involves interviewing each patient and performing necessary tests on a case-by-case basis.In the KURAMA cohort, all physicians used a common medical questionnaire to record patient states over time under certain standards and quality controls.It is necessary to standardize data and improve the medical system to enable the collection of high-quality time-series data that are useful for both medical care and analysis while balancing routine medical care and data analysis and eliminating unnecessary burdens as much as possible [58,59].

Conclusions
This study suggested that evaluating state multistability and determining the patient's state in daily practice may enable treatment plan optimization over the entire course of treatment.We believe that this study will contribute to the development of personalized medicine utilizing real-world data.

Supporting information
S1 Checklist.STROBE statement-checklist of items that should be included in reports of observational studies.(DOC) After the inspection is completed, the seven factors are entered into a tablet by the physician, and the state number, energy, and coordinate of the multistability score in the state transition are displayed.In this example, the energy is higher than the threshold value and is in a transitional state; therefore, the prospects for remission are considered adequate.However, the patient had a poor-stability cluster and is currently using csDMARDs.A change to a more effective drug may be considered.https://doi.org/10.1371/journal.pone.0302308.g007

Fig 1 .
Fig 1. Patient state transition model in RA drug treatment.The higher the energy is, the more unstable the state and the easier it is to obtain a therapeutic effect.At initial consultation, many patients present in a poor state and with high disease activity (i.e., high energy).With treatment, they may transition to a good state or stabilize in a poor state.The physician observes the patient's state and executes a treatment strategy to stabilize the patient in a good state or to prevent stabilization in a poor state below a certain energy threshold.The objective of this study was to visualize the state transitions of these patients as a population.
Fig 4C shows the threshold value of the energy at which pattern transitions can occur (see: "Disconnectivity graph").

Fig
Fig 4D shows the distributions of state numbers and energies.The horizontal axis is the energy value, and the vertical axis is the number of states.The line in red indicates the energy threshold (-1.48) at which the pattern can move between states.The threshold value was approximately -1.48.Approximately 85% (109/128) of the state numbers had energy values that could be transferred between patterns.The state numbers were aggregated into four quadrants at the energy threshold, and the number of patients per quadrant was calculated (Fig 5).Fig 5A divides energy and pattern into four quadrants by threshold.The quadrant in which the energy is lower than the threshold and the activity pattern is "good stability" is G-L; the quadrant in which the energy is higher than the threshold and the activity pattern is "good stability" is G-H; the quadrant in which the energy is higher than the threshold and the activity pattern is "poor stability" is P-H; and the quadrant in which the energy is lower than the threshold and the activity pattern is "poor stability" is P-L.Fig5Bshows the number of patients in each quadrant, with a marked increase in G-L and a decrease in P-L until approximately two points (i.e., 6 months), followed by a gradual increase or decrease.P-L decreased until approximately two points (i.e., 6 months), with no significant change thereafter.G-H decreased significantly throughout the entire period.G-H did not increase or decrease significantly across the entire treatment period.Fig5Cshows the HAQ scores in each quadrant: G-L patients met the remission criteria for almost the entire period; the number of P-L patients who met the remission criteria decreased over time; the number of G-H and P-H patients who met the remission criteria changed

Fig 7 .
Fig 7.Examples of future real-world clinical implementations.After the inspection is completed, the seven factors are entered into a tablet by the physician, and the state number, energy, and coordinate of the multistability score in the state transition are displayed.In this example, the energy is higher than the threshold value and is in a transitional state; therefore, the prospects for remission are considered adequate.However, the patient had a poor-stability cluster and is currently using csDMARDs.A change to a more effective drug may be considered.